# analiza kowariancji
library(faraway)


data(sexab)
sexab
help(sexab)

# zwykle roznice miedzy grupami i prezentacja
by(sexab,sexab$csa,summary) 
plot(ptsd~csa,sexab)
plot(ptsd~cpa,pch=as.character(csa),sexab)
library(lattice)
xyplot(ptsd~cpa|csa,sexab)
xyplot(ptsd~cpa|csa,sexab,type=c("p","r"))

t.test(sexab$ptsd[1:45],sexab$ptsd[46:76]) # istotna stat.


g=lm(ptsd~cpa+csa+cpa:csa,sexab)
summary(g)  # brak istotnosci roznicy w "nachyleniu" pr. regresji

sexab$csa=relevel(sexab$csa,ref="NotAbused")  # zmiana odniesienia
g=lm(ptsd~cpa+csa,sexab)

summary(g)


plot(ptsd~cpa,pch=as.character(csa),sexab)
abline(coef(g)[1],coef(g)[2])
abline(coef(g)[1]+coef(g)[3],coef(g)[2],lty=2)


plot(g,1:6)
plot(fitted(g),residuals(g),pch=as.character(sexab$csa))


###############  przykład z wieloma poziomami


data(fruitfly)
help(fruitfly)
summary(fruitfly)
plot(longevity~thorax,fruitfly,pch=unclass(activity))
legend(0.63,100,levels(fruitfly$activity),pch=1:5)

xyplot(longevity~thorax|activity,fruitfly,type=c("r","p"))


g=lm(longevity~thorax*activity,fruitfly)
summary(g)
plot(g)
anova(g) # F test
step(g) #ten sam rezultat
gb=lm(longevity~thorax+activity,fruitfly)
summary(gb)

plot(residuals(gb)~fitted(gb),pch=unclass(fruitfly$activity))
plot(gb)

#xyplot(residuals(gb)~fitted(gb)|activity,fruitfly)

library(MASS)
boxcox(gb)
gt=lm(log(longevity)~thorax+activity,fruitfly)
summary(gt)
exp(-0.12) # tylokrotnie skraca się czas życia
exp(-0.42)
exp(2.72146)^0.1
plot(residuals(gt)~fitted(gt),pch=unclass(fruitfly$activity))
plot(gt)


#######################################
#### jednoczynnikowa analiza wariancji


data(coagulation)
summary(coagulation)
help(coagulation)
g=lm(coag~diet,coagulation) # standardowo w R przyjmujemy alpha_0=0;
summary(g)  #p-value mowi o istotnosci sredniej w dieta A, a sredniej w obecnej diecie
gi=lm(coag~diet-1,coagulation) # teraz \mu =0, a alpha_i oznaczaja srednie w kolenjnych grupach (A,B,C,D);
summary(gi) #uwaga - R^2 dla takich modeli jest źle obliczony
gnull=lm(coag~1,coagulation)
#mean(coagulation$coag)

anova(gi,gnull)  # porownanie istotnosci takiego modelu, a uzycia sredniej
#anova(g,gnull) to samo
plot(g)
# test Levene'a
med=with(coagulation,tapply(coag,diet,median))
med
ar=with(coagulation,abs(coag-med[diet]))
ar
anova(lm(ar~diet,coagulation))


TukeyHSD(aov(coag~diet,coagulation))




model<-aov(coag~diet,coagulation)
summary(model)  ### the difference between means is sufficient

plot(model)

library(stats)
library(car)
#var.test(CD14.2191~Mutation,AML,ratio=1)
fligner.test(coag~diet,coagulation)  #Flignera-Killeena
bartlett.test(coag~diet,coagulation) #Bartletta, leszy dla normalnych
levene.test(coag~diet,coagulation)
leveneTest(coag~diet,coagulation)  # dosc odporny na normalnosc



TukeyHSD(model) 
library(agricolae)
HSD.test(model,"diet",group=TRUE,console=TRUE) ## the number of Male and Female should be equal!!!
SNK.test(model,"diet",group=TRUE) ## the number of Male and Female should be equal!!!
LSD.test(model,"diet",group=TRUE)
scheffe.test(model,"diet",group=TRUE)


boxplot(coag~diet,coagulation)

#############################
library(PBImisc)
#AML$change=AML$CD14.2191-AML$CD14.control
#attach(AML)
#AML
#by(change,Mutation,mean)
 
#boxplot(change~Mutation,AML)
help(AML)
summary(AML)
summary(aov(CD14.2191~Mutation,AML))

model<-aov(CD14.2191~Mutation,AML)
summary(model)  ### the difference between means is sufficient
plot(model)

library(stats)
#var.test(CD14.2191~Mutation,AML,ratio=1)
fligner.test(CD14.2191~Mutation,AML)  #Flignera-Killeena
bartlett.test(CD14.2191~Mutation,AML) #Bartletta, leszy dla normalnych
levene.test(model,center=median)
leveneTest(model,center=median)  # dosc odporny na normalnosc



TukeyHSD(model) ## the number of Male and Female should be equal!!!
library(agricolae)
HSD.test(model,"Mutation",group=TRUE) ## the number of Male and Female should be equal!!!
SNK.test(model,"Mutation",group=TRUE) ## the number of Male and Female should be equal!!!
LSD.test(model,"Mutation",group=TRUE)
scheffe.test(model,"Mutation",group=TRUE)